Simulation and Correction of the Beam 
Hardening effect in X-ray Tomography 


by 

K. Rama Krishna 


TH 

k857^ 


NUCLEAR ENGINEERING AND TECHNOLOGY PROGRAMME 

INDIAN INSTITUTE OF TECHNOLOGY KANPUR 

Aoril. 2000 




Simulation and Correction of the Beam Hardening effect 

in X- ray Tomography 


A Thesis Submitted 

in Partial Fulfillment of the Requirements 
for the Degree of 

Master of Technology 


by 

K. Rama Krishna 



Nuclear Engineering and 1 echnology Programme 
Indian Institute of Technology, Kanpur 
April 2000 



Certificate 



It is certified that work contained in the thesis entitled. Simulation and correction of 
the Beam Hardening effect in X - ray tomography, by K. Rama Krishna, has been 
carried out under our supervision and this work has not been submitted elsewhere for 
a degree. 


Dr. Prabhat Mun.shi, 

Professor, 

Department of Mechanical 
Pngincering, 

Indian Institute of 'feclinology, 
Kanpur 208 016. 


Dr. M. S. Kalra, 

Professor, 

Department of Mechanical 
Engineering, 

Indian Institute of Technology, 
Kanpur 208 016. 


April, 2000. 



Acknowledgements 


, ^ r u ■ pratitude and appreciation to my 

I would like to express my deep felt sincere grau 

. T\r»x/r i- jTxx^c tCalra for their skillful guidance, 
thesis supervisors Dr. P. Munshi and Dr. M. o- 

constant supervision, timely suggestions in carrying ont the present work 

, • , , ^ oHHhar for his tireless supervision, 

I owe my sincere thanks to Dr. K. Muralicinar 

encouragement and moral support throughout the course of this work. 

to my beloved parents for 

I acknowledge my heartiest regards and reveren 
their deep affection and moral support. 

t r ■ j c * n'prhan Sharad and Ruchi for helping 

I would like to thank my friends Suneet, l-ariian, 

j -I- 1, 1 ♦ u TiiP'ir company made my stay at iitk 

me and providing a homely atmosphere, iheir cun i 

memorable. 



Abstract 


Polychromatic X-ray sources are universally used in computerized 
tomography to obtain adequate intensity. These sources, however, can produce some 
artifacts in the reconstructed image due to non-linearities. Beam hardening is one 
such artifact, which produces false line integrals. 

I’he present investigation involves how one can estimate from the total 
attenuation, p, of a polyenergetic X-ray beam what the total attenuation, m, of a 
monoenergctic beam would have been along the same ray. In the present work three 
different specimens are studied which have different geometries. Simulation results 
show the introduction of beam hardening although the projection data are 
qualitatively same but quantitatively different. 

The polynomial approximation method is working well. For a typical X-ray 
spectra passing through the material one can find a simple function / such that f(p) 
is a sufficiently close estimate of m to allow good reconstructions. 
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Chapter 1 


Introduction 


The technique of computerized tomography(CT) has established itself as a leading 
tool in diagnostic radiology over the past thirty years and is catching on fast in the 
nondestructive evaluation area in a variety of situations. The first CT scanner was 
developed by Hounsfield (1973) who was later on awarded the Nobel prize for medicine. 
In past, CT was done by X-rays but now other sources like y-ray, laser, ultrasonic 
resonance etc. arc also used, 'fhere has been a very great deal of activity in recent years to 
find algorithms which arc fast when implemented on a computer and which produce 
acceptable recoirstructions in spite of the finite and inaccurate nature of the data. 

1.1 An Overview of computerized tomography 


'fhe aim of computerized tomography is to assign to every point inside the object 
a number, which is specific to the material occupying that point. A suitable candidate for 
this number is the X- ray attenuation coefficient of the material. 



‘oui 


SOURCB 


OBJECT 

f(x,y) 


DETECTOR 


Mg. 1 . 1 Different materials attenuate the X- ray differently. 



Inlroduclion 


The object function is denoted by fix,y). The relation between the intensities 
7,^ and in Fig. 1.1 is described by 


L 

By taking the natural logarithm we get 


( 1 . 1 ) 


in 



out 




( 1 . 2 ) 


This means that what we read in detector is the line integral of the object function. 
This object function f{x,y), for X-ray / y -ray setup, represents the linear attenuation 
coefficient, //(.v. v), characteristic of the material. 



Fig. 1 .2 A Radon value is the line integral along the line L. i 

The radon transform R[f (r, <?)], or the projection p(r, ^)j, of an object /(x, y) | 

is the line integrals through the object in all possible directions. This means, a single 
Radon value is the integral of all points along a line with angle and perpendicular 
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distance r from the origin (Fig 1.2). The Radon transform (Herman 1980) can be written 
as 


p{r,0) = R [/(r,<9)]= J|/(x,y) <5(xcos(9 + j^sin6'-r) dx dy (1.3) 
An important relationship between the Radon transform and the Fourier transform 
is the Fourier slice theorem. The theorem states that the one dimensional Fourier 
transform of the Radon transform along a radial line is identical to the same line in the 



Fig. 1 .3 The Fourier slice theorem. 
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two dimensional Fourier transform of the object. 

F [R [/(i?,^)]] = [/(i? cos Sint?)]] (1.4) 

The theorem in illustrated in Fig. 1.3 . 



_ , Reconstruction 

Radon space 


Fig. 1 .4 Computing the inverse of the Radon transform. 

1.2 Image Reconstruction from projections 

'fhe range of applicability of image reconstruction is very wide. At one end, data 
from electron microscope are used to reconstruct the molecular structure of 
bacteriophages, while at the other end data collected by the rockets sent outside the 
earth’s atmosphere arc used tp reconstruct the X- ray structure of the supernova remnants. 
Image reconstruction can very aptly defined in words of Herman (1980) as : 

“ Image reconstruction from projections is the process of producing an image of a 
two dimensional distribution (usually of some physical property) from estimates of its 
line integrals along a finite number of lines of known locations”. 
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1.3 An overview of the present work 

The measurements in computed tomography can only be used to estimate the line 
integrals. Inaccuracies in these estimates are due to the width of the X-ray beam, 
hardening of the beam, photon statistics, etc. Radon’s inverse formula is sensitive to these 
inaccuracies. 

When an X-ray beanr passes through the material, its attenuation at any point 
depends on the material at that point and on the energy distribution (spectrum) of the 
beam. A difficulty arises due to the fact that the X-ray beam used in computed 
tomogarphy consists of photons at different energies. The attenuation at a fixed point is 
generally greater for photons of lower energy as a result the energy distribution 
(spectrum) of the X-ray beam changes (hardens) as it passes through the material. X-ray 
beam reaching a particular point inside the material from different directions are likely to 
have different spectra (having passed through different materials before reaching the 
point of interest) and thus will be attenuated differently at that point. This makes it 
difficult to assign a single value for the attenuation coefficient at that point inside the 
material. 


1.3.1 The problem of Beam Hardening 

I'hc aim of computed tomography is to obtain information regarding the nature of 
material occupying exact positions inside the object. 

In computed tomography we have two sets of measurements: 

(i) Calibration measurements, on which we can base an estimate of what the detector 
measurement.s would be if the object to be reconstructed is not between the source and 


the detector . 
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(ii) Actual detector measurements with the object of interest in position. 

The method by which these measurements are taken is shown in Fig. 1.5. There is a 
region of space (referred to as reconstruction region in Fig. 1.5 ) is occupied by some 
homogeneous reference material such as air or water during the calibration measurement. 
During the actual measurement, the object of interest is inserted into the reconstruction 
region, (partially) replacing the reference material. It is an important restriction that the 
object of interest does not occupy any point outside the reconstruction region. 

Suppose that we have a monoenergetic X-ray source with photon energy e . For 

a 1 ixed source and detector pair, let be the calibration measurement, and let be 
the actual measurement. We can define monochromatic ray sum, m, for this beam by 


m 


- In 


'A' 


(1.5) 


and we refer to the set of m's for all source and detector pair positions as the 
monochromatic projection data. In practice, the X- ray beam is polychromatic. Let 

and denote calibration and the actual measurement, respectively. We can define 
polychromatic ray sum, /?, for this beam by 


P 


- In 


J 


( 1 . 6 ) 


and we refer to the set p’s for all sourcc-dctcctor pair positions as the polychromatic 
projection data. 

The problem of beam hardening is that for any source and detector pair we can 
obtain p, but reconstruction procedure requires m as per Eq. (1.5). Beam hardening 
re.sull.s in false gradients of the linear attenuation coefficients in the CT cross section 
images, indicating a non-existent density or composition gradient in the imaged object. 
Convciion for beam hardening effect is a must in Quantitative Computerized 
Tomography (QCT). For a given total attenuation of polyenergetic X- ray beams 
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through the object, one has to estimate the total attenuation of monoenergetic X — ray 
beams through the same parts of the object which are precise enough for useful 
reconstruction of the monoenergetic linear attenuation coefficients in the material object. 
The mathematical and computational procedures in CT imaging are summarized in Fig. 
1.5. 



Fig. 1.5 


Outline of the mathematical and computational procedures in Cl'. 
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1.4 Thesis layout 


• Chapter 2 gives the detail of the reconstruction algorithm (CBP) and 
derivation of formulae for polyenergetic radiation. 

• Chapter 3 gives the details of the simulation beam hardening effect of the 
specimens. 

• Chapter 4 gives the details of the method for removing the beam hardening 
effect in computerized tomography. 

• Chapter 5 provides discussions on the results based on the simulated data. 

• Chapter 6 summarizes the conclusions and suggests scope for further work. 


8 
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Chapter 2 


Theoretical Formulation 


In this section, the computational and mathematical procedures underlying the 
data collection, image reconstruction, formula for polyenergetic radiation and image 
display used in the C'F are discussed. 


2.1 Preliminaries 

TIic number of counts of photon after passing tlirough a curve ‘c’ in the material 
being te.stcd is given by 

N = NoQ,x^(r\^^{r,(f))dl)) (2.1) 

where the integration is along the chord length of ‘c’, 

N is the number of photon counts after traversing the chord length, 
i.s the initial photon counts and 
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jj, is the attenuation coefficient. Since ^ depends on the material and the energy of the 
radiation, a distribution of ^ is a direct indicator of the material composition of the 
component being studied. 

Equation 2. 1 can be written as 


In 






= p{s,e) 


( 2 . 2 ) 


where pi^s^O) is called the projection data for the tomographic algorithm, and it is the 
integral of the function along the line specified by S and 0 (Fig. 2.1). 

The aim of tomograph}- is to reconstruct the function ^) , if a set of several p-values 

(p{s\0)) measured along various chords (c) is given. This is the fundamental problem 
of Cl' and CBP has been used in the present study for that purpose. The p values can be 
suitably normalized to get the material density distribution, if so desired. 


2.2 Data Collection Mode 


'fhc image proccs.sing methodology requires attenuation data to be collected by an 
array of radiation detectors for the reconstruction of the function pir^tj))- In this study 

the mode of collection is the Parallel Beam Geometry (PBG) mode (Fig. 2.1). This 
system consists of several pairs of radiation source and radiation detector, which can scan 
the object completely. The SD (source detector) pairs are spaced uniformly and the object 
can be rotated to get the data for different views. The line SD represents the path of the 
data ray or the chord along which the function p{s^6) can be found out. The 


in 
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perpendicular distance from the center of the object to the path of the ray is denoted by 
S . The object table is rotated to get several sets of 'p' for different values of 6 . 



Fig. 2.1 Data collection geometry. 

Typical values of the number of rays and the number of projections used in the 
present work arc 100 and 100. The source-detector system is moved pixel by pixel, thus 
generating a grid of 100 by 100 over the physical region of interest. 


LfBRAlP 

I I T., nummn 


r. A 130122 
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2.3 CBP Algorithm 

Convolution Backproj action (CBP) algorithm has been described in great detail 
by Herman (1980). In this section we review CBP briefly as reported earlier by Munshi 
(1992). 

Fig. 2. 1 shows the data collection geometry for a parallel beam CT scanner. The 
object function, denoted for generality by in subsequent equations, is 

represented by a unit circle and one (of many) data rays is represented by SD. The ray 
indices are S and 0 , where S is the perpendicular distance of the ray from the object 
center, and 0 is the angle of the source position (or object rotation). The CT data denoted 
by (?) given by 

p(s,e)= y^fir,,/))dz (2.3) 

1 lore, Z is the variable of integration along the chord SD. The CT machine collects the 
projection data p(^s,0} for many values of ^ and 0 . The “projection slice theorem” 

Herman(1980) states the equivalence of the two-dimensional Fourier transform of 
f{r,(^) and the 1 -dimensional Fourier transform of pis, 6) with respect to S. 

Symbolically 


p{R, 0) = /(?? cos e, R sin 6) (2.4) 

where the symbol A represents the Fourier transform and R is the spatial Fourier 
frequency. A two-dimensional Fourier inversion of Equation 2.4 leads to the well-known 
tomographic inversion formula 

/(^(^) = ] i ^ I (2-5) 

O-oo 


where 
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p{R,e) = ^^p{s,0)e-‘^’^ds (2.6) 


We note that the inner integral in Equation 2.5 is divergent. A practical implementation 
of Equation 2.5 incorporates the replacement of the factor | i? | by | i? | fF {R) . Where 
JV (R) is a suitable window function that vanishes outside the interval Here, 

I I is the Fourier cutoff frequency. Normally, W{R) is an even function of R . Thus 
Equation 2.5 takes the approximate form 
TT 00 

f{r,(f)) « j I R \ W{R)dRde (2.7) 

0— 00 

If p(R^0) also vanishes for | 7? [> 7?^, then the reconstructed function , denoted by f, 
agrees exactly with /'(r, tp) with the following window function 

= 0, I I> .R. (2.8) 

Equation 2.8 and the convolution theorem of Fourier transforms give the reconstructed 
function f as 


TT 00 

/('■. <t) = \ \ P.9) 

O-oo 

where 


and 


00 


q(s)= l\R\W{R)e‘^’^^dR 


-00 


s' = r cos(0 - <p) . 


( 2 . 10 ) 


( 2 . 11 ) 


n 
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The index s' is the data ray passing through (r,^), the point being 
reconstructed. The inner integral in the Equation 2.9 is a one-dimensional convolution 
and the outer integral, corresponding to the averaging operation (over 0 ), is termed as 
backprojection and hence the name convolution backprojection for this particular 
implementation. The CBP method is also known as the filtered backprojection algorithm 
because of the “filtering” of the Fourier transform of the projection data p by the 

window (or filter) W (/?) in the initial stages of the formulation given by Equation 2.5. 

The function Cj[{s) , known as the convolving function, is evaluated once and stored for 

the repeated use for different views (or different angles 6 ). 

For a given point (r, ^), the inherent error in the CBP implementation. 
Equation 2.9 is, 

fir, (()) ( 2 . 12 ) 

where f and J are given by Equation 2.5 and Equation 2.9 respectively. This error is 
strictly due to finite cut-off, of the Fourier frequency and is precisely zero if the 
projection data happens to the band-limited and the cut-off frequency is chosen to be 
cho.sen to be the highest frequency contained in In general, to avoid aliasing artifacts 

Merman ( 1 980) has recommended the choice, = l/(2A.y), where Iss is the spacing of 
the data rays. 

2.4 Formula for polyenergetic radiation 

'fhe linear X-ray attenuation coefficient at a point inside a cross section of the 
object depends on the position of the point (x,y) and on the X- ray energy e. It can be ; 
denoted as //(XjT'td). From equation 2.2 in case of monochromatic beam it can be 

! 

written as ! 


14 
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( 2 . 12 ) 


In case of polychromatic beam the result will not be , but rather an estimate for the 
more complicated integral 


p, = -In Jr(e) exp - jp(x,y,e) 

0 K L 


de . 


(2.13) 


Where r(t’) is the probability that the detected photon is at energy e . 

It is derived in the following mamter (Herman (1979)): 

Let denote the linear attenuation coefficient at energy € at the point Z 

between the source {Z Zs< -■= 0) and the detector (Z = Zd> = D). 
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Let Pg {z) is the probability that a photon of energy € which leaves the source in 
the direction of the detector will get as for as the point Z without being removed from the 
beam. Let (z, Az) is the probability that a photon of energy e, which already traveled 

up to Z from source is removed before it gets to a distance Z+AZ from the source. Then 
we can write 

p,{z + Az) = pfz) [l-e,(^,Az)] 

We can write (z) as 


(2.14) 


Az 


(2.15) 


From Hq. 2.14 


p Pz + Az)-p^,iz) 1 q,iz,^z) 

p, (z) Az 

Taking A/,---)-!), wc get 

Pe (^) ^ 


From liq.s. 2.15 and 2.17 


P[.(z) 

pA^) 


/L-(2)- 


(2.16) 


(2.17) 


(2.19) 


Integration of both sides from Z, to Z,, leads to 

In p, iz,i ) - In p, (z, ) = - \m, {z)dz . (2.20) 

Recalling the definition of /7,(z), = 0), and so 
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*c/ 

= - \jufz)dz. 


( 2 . 21 ) 


pfzj) = exp 



( 2 . 22 ) 


Suppose that S^. photons of energy e leave the source in the direction of the 
detector, and that 6). is the efficiency of the detector at energy 6 . Then, on an average, 
the number of photons of energy e counted by the detector is 




where 


exp 


\p,{z)dz 


= exp 


1 1 

\p,{^)dz 


L 0 


(2.23) 


= exp 


w 

\Mfz)dz 


exp 


\p,{z)dz 


(2.24) 


Ily the Hq. 2.24 the number of photons at energy B counted by the detector during the 
calibration mea.surement is 

/> 


= S^S^k^ exp 


\p°{z)dz 


(2.25) 


Wlierc is the linear attenuation coefficient of the reference material. Let is the 

total number of photons detected during the calibration measurement, then the detected 
spectrum is 



(2.26) 


llencc polychromatic ray sum is 


1 ^ 
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p = - In 




\^p J 


-In 


a, 

js,SX exp 


~\pXz)dz 


de 


■In 


cr, exp 




c. 


de 


In Jr..(--)cxp-j(//,(z)-/.:(z) )dz 


L 0 


de 


(2.27) 


If reference material is air then 0. 


2.5 Displaying of an image 


Fnr displaying of an image the program written by Farhan (1999) is used. For the 
purpose of displaying the image, the CT numbers are read from the CBP output file and 
apj’uupfiatc grey levels a.ssigned corresponding to it. Thus by generating pixels for each 
of the projection data at their position the CT image can be graphically displayed. The 
color code is printed adjacent to it for a quick reference. 

I'hc color code is a must to extract information from the colored images produced 
in the present w'ork. 'I'hc natural choice for the color code was our own solar spectrum 
“VlBtiYOR” in the scale of 0 to 255. Number 0 is represented by ‘Violef and 255 by 
‘Red’ and the intermediate colors linearly interpolated to represent other numbers. 
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2.5.1 VIBGYOR Representation of CT Images 

The colours in the image are assigned in a manner such that each distinct value in 
the reconstructed matrix is displayed in a distinct colour. The colour of the pixel is 
determined by creating a linear scale such that the smallest number in the matrix gets the 
colour violet and the largest number gets the colour red. The intermediate colours are 
generated by linearly interpolating between the RGB (Red, Green, and Blue) values of 
the colours of a rainbow (VIBGYOR: Violet, Indigo, Blue, Green, Yellow, Orange, and 
Red). 'I’he program is written with tlie OpenGL library in UNIX platform. Reconstruction 
of tomographic data results in a matrix in which each element corresponds to the 
magnitude, of the material property being observed, at that point. If we are able to view 
this ntatrix in the form of a coloured image, each colour in the image will correspond to 
the different observed values of the material property, thus providing us with a view to 
the interior of the object under observation. The VIBGYOR pattern provides a good 
mode of \isuali 7 iug the variation in the material properties of the object. This 
\ isuaIi/ation had been extremely simple if we had a fixed number of values for the 
ob.scivable material property and if we had a unique identifier for each colour in a 
VIBGYOR spectrum. Thus the problem can be broken up into two pieces: 

1. Identifying the RGB values for each of the colours required along the VIBGYOR 
spectrum. 

2. .Xssi-Miini' tiie elements of the matrix to one of these colours based on a linear scaling 
of all the values in the matrix. 

It is a common practice to represent a colour, as a set of three numbers, which 
coijc-qHMh! to tlic intensities of Red, Green and Blue, respectively, that when mixed 
would prtiduce the colour. To tackle the first sub-problem, we took the RGB values of the 
seven main colours of the rainbow as: 

« 


Violet : 148, 0,211 



Theoretical formulation 


• Indigo ; 75, 0, 130 

• Blue : 0, 0, 255 

• Green ; 0, 255, 0 

• Yellow: 255,255, 0 

• Orange; 255, 165, 0 

• Red : 255, 0, 0 

In each case 0 corresponds to the minimum possible intensity and 255 
corresponds to the maximum possible intensity of the colour (Red, Green or Blue). The 
program prompts the user to enter the number of colours desired along the VIBGYOR 
spectrum. Keeping the values of the major colours as given above, the other values are 
ctiinputed by linearly interpolating the Red, Green and Blue components between these 
major colours. 11" 20 colours are required to lie between Blue and Green, then the Red 
component of till these colours is assigned a value of zero. The Green component is 
calculated as Cr>ce«,- Ti ( ((2J5-(9) x / /20), and similarly the Blue component of these 
colours is calculated as Bluer 25 5 +((0-255) x i/ 20). 

The second sub-problem is also tackled by using linear interpolation. The 
minimum and maximum values occurring in the matrix are determined first. If the 
number ol'coUmrs required by the user is n, then n-2 values lying evenly spaced between 
the two limiting values arc computed. Now we have n values corresponding to each of 
the n required colours. While di, splaying the image the colour to be assigned to each pixel 
i.H determined by comparison with these n values. This can be better understood with the 
following example. Suppose the user asks for a 256 colour long VIBGYOR spectrum, let 
the smallCsSt value found in the matrix be Vmin and the largest value is V^ax- An array of 
256 numbers i.s created and Vmin and Vmax are placed in its first and last place, 
respectively. These and the other 254 values can be computed using the formula; 


VrV„,n + {V^,rVmm)^ i/256 


0<i<255. 
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Chapter 3 


Simulation of beam hardening effect 


in thi.s chapter simulation procedure of beam hardening effect is described for 
assumed cross sections with the help of chosen X-ray spectra. 


3,1 Simulation of polyenergetic projection data 


'Fhe formula used for simulation of polyenergetic projection data is derived (Eq. 
2.2K) in chapter 2, It is assumed that the spectrum of the X-ray beam can be 

approximated by a discrete spectrum consisting of J energies e(l), e(2) , e{J) and 

that is the probability that a detected photon(of the X- ray beam through air 

bctivcen the source and the detector ) is at energy e(J) . Here it is assumed that air is the 
only reference medium so, that /r" = 0 and the equation become simpler. Let us divide 

the cro.ss section into / pixels. We will try to estimate the linear attenuation coefficient in 
each of the / pixels. From the equations 2. 12 and 2.27 of the chapter 2 we have 
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m = 




(3.1) 


P = 



(3.2) 


where Z' denote the length of intersection with the zth pixel of the line from the center of 
the source to the center of the detector and is the linear attenuation coefficient at 

energy v in the zth pixel. Simulation is done on three different types of cross sections. 
I’irst specimen (SI) is a cross section consisting of three small circles maintained at the 
same density at each of the energies. For this specimen the X-ray spectrum is taken from 
llennan(197')). 


.i 


circles f,i 
(in cm"’) 

non circles p, 

(in cm"') 

Energy 

(keV) 

1 

O.I 

0.999 

0.265 

41 

2 

0.4 

0.595 

0.226 

52 

.3 

0.3 

0.416 

0.210 

60 

4 

0.2 

0.208 

0.174 

100 



Table 3 . 1 X- ray spectrum for SI. 





Simulation of beam hardening effect 



Min = 0.0000 Max =0.4160 

LAvg= 0.2512 AAvg = 0.2404 

File N Ray : 100, NView : 100 


Fig. 3.1 Original image of SI at energy 60 keV. 



Min = 0.0579 Max = 0.4961 
LAvg = 0.2666 AAvg = 0.2545 
File : ,NRay : 100, NView : 100 


Fig. 3.2 Beam hardening reconstructed image of S 1 . 


Simulation of beam hardening effect 

In original image (Fig. 3.1) //= 0.210 is represented as green colour and 
fx = 0.416 is represented as red colour. After beam hardening effect (Fig. 3.2) maximum 
/d found out to be 0.4961. The images are normalized between minimum 0 and 
maximum 0.4961. 

Ihc second specimen (S2) is a cross section consisting of two circles, two 
triangles, two rectangles. The materials inside them are iron, titanium and void. The 
spectrum is taken from a NIST (National Institute of Standards and Technologies) 
I lubbelK 1 9H2). 'flic density of void is taken as zero at all the energies. 


.i 


Iron }.i 
(in cm"’) 

Titanium |u 
(in cm"') 

Energy 

(l*10"’MeV) 

I 

0..3 

2.926 

1.235 

1 

2 

0.4 

1.1496 

0.596 

2 

3 

0.3 

0.8653 

0.473 

3 


Table .3.2 X -ray spectrum for S2. 


Fig. .3.3 is original image of S2 in which titanium is shown in blue colour 
(/r 0.473 ), iron is in yellow colour (/r = 0.8653) and void is in violet colour. After the 

introductiim of the beam hardening effect the maximum pixel value becomes 1.3191. The 
images are iUTmali/cd kdween minimum 0 and maximum 1.3191. 



Simulation of beam hardening effect 



Min = 0.0000 Max = 0.8653 
LAvg = 0.5199 AAvg= 0.4962 
File ;,NRay : 100, NView : 100 


P'ig. 3.3 Original image of S2 at energy 0.3 MeV. 



Min = 0.0201 Max = 1.3191 
LAvg = 0.7463 AAvg = 0.7099 
File : ,NRay : 100, NView : 100 


Fig. 3.4 Beam hardening reconstructed image of S2. 


Simulation of beam hardening effect 


Third specimen (S3) is a star type of object whose /u values at different energies 
are shown in Table 3.3 which is taken from Herman (1979). There are eight small circles 
representing the void whose fj. values are taken as zero at all the energies. 


j 


star \x 

(in cm"’) 

non star area \x 

(in cm”') 

Energy 

(keV) 

1 

0.4 

0.999 

0.265 

41 


0.3 

0.595 

0.226 

52 

3 

0.1 

0.416 

f 

0.210 

60 

4 

0.1 

0.265 

0.183 

84 

5 

0.1 

0.208 

0.174 

100 


Table 3.3 X-ray spectrum for S3. 

In the original image (I’ig. 3.5) the maximum pixel value is 0.416. After the beam 
hardening effect (Fig. 3.6) the maximum pixel value changed to 0.6775. The linear 
attenuation coefficients of voids arc also changed. The images are normalized between 
minimum 0 and nuixlmum 0.6775. 



Simulation of beam hardening effect 



Fig. 3.5 Original image of S3 at energy 60 keV. 



Min = 0.0925 Max = 0.6775 
LAvg = 0.3397 AAvg = 0.4922 
Fite :,NRay : 100, NView : 100 


Fig. 3.6 Beam hardening reconstructed image of S3. 
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Chapter 4 


Correction of beam hardening effect 


4.1 Correction method for beam hardening 


In this chapter the correction procedure is summarized. For correcting the beam 
hardening effect the procedure adopted by Herman (1979) with few a modifications is 
used. 

For a fixed source and detector position, let m denote the monochromatic ray 
sum (the measured integral of the linear attenuation coefficient between source and 
detector if we had an X~ray source that produces only X-rays at energy e ) and let p 
denote the polychromatic ray sum (estimate of the integral of the linear attenuation 



Correction of beam hardening effect 

coefficient between the source and the detector based on measurements with actual 
polychromatic X— ray source). Using our earlier notation m and p can be defined by 

( 4 . 1 ) 


m = ^/d{x,y,e)dl 


p = -\n Jz-(e) exp - ^p{x,y,e) de . (4.2) 

0 K i j 

Here E is the highest energy level present in the beam, and r(e) is the probability that a 
detected photon (of the X-ray beam through air between source and detector) is at energy 
<?. 


'flic least expensive type of beam hardening correction can be done by using a 
function / which is such that, for each source/detector pair, f(p) is a reasonable 
estimate of m. Let us refer to the reconstruction from the so corrected polychromatic 
data { f(p) } as the first reconstruction. It is a set of I number, p[ , representing the 

estimated linear attenuation coefficient at energy e of the material in the zth of a total 
of /pixels. 

It is assumed that the spectrum of the X-ray can be approximated by a discrete 
spectrum consisting of /energies e(l),e(2) ....e(/) and that is the probability 
that a detected photon is at energy e{j ) . We further assume that there are functions 
, ^ 2 ,. . . gj such that, for !</</, g^ ) is a reasonably good estimate of the linear 
attenuation coefficient at energy e{j) of the material in the ith pixel. We describe one 
such set of functions below. We use to denote g^ {p[ ) . For a fixed source detector 

pair, let Z' denote the length of intersection with the ith pixel of the line from the center 
of the source to the center of the detector. We define monochromatic pseudo ray sum 
fn , and polychromatic pseudo ray sum p by 
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m = 


I/-.; Z' 



J 


J = 1 




(4.3) 


(4.4) 


We see that m approximates to m, and p approximates to p, and hence f (p) 

approximates iof{p). Furthermore, since the line integrals in Eqs. 4.1 and 4.2 are 
approximated in the same way in Eqs. 4.3 and 4.4, it appears likely that the errors 
m — m and f {p) — f {p) in these approximations will be similar (i.e., the difference 
between these errors will be considerably smaller than either of the errors ). Hence 
~ J (p) + f (p) is an approximation to m is superior to the use of just f (/?) . This 
is true in the sense that 

A({/(/7) + m- f{p)}, {m}) < A({/(p)}, {m}) (4.5) 

where A represents the root mean square error. The second reconstruction is one 
obtained from the data w -/(p)4- f{p). Since the second reconstruction is 
presumably more accurate than the first one, this process can be repeated 

The correction procedure followed in the present work is shown Fig. 4.1. For 
polynomial curve fitting Matlab is used. Matlab solves the system of linear equations in 
order to obtain the coefficients of the polynomial. The projection data p, is a matrix 
comprising of rows representing number of views and columns representing number of 
rays. The curve-fitting formula relating p and m is 

[m]^c^+c^[p] + c^[pf+ + c^[pY (4.6) 

Where the coefficients C- ' $ arc found by the generalized vandermonde matrix method. 



Correction of beam hardening effect 
For the simple case 


y = CQ +c^x + C2X^ + + c^x" 

this method reduces to solving the system of linear equations as 



(4.7) 


(4.8) 


For regression analysis using matrices, this is expanded into a matroid equatiqn. 

For the specimen S 1 the initial guess is taken as cross section consisting of two 
circles. For the specimen S2 the initial guess is taken as cross section consisting of one 
triangle and one circle. For the specimen S3 the initial guess is taken as cross section 
consisting of single circle at the center. 
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Fig. 4. 1 Flow chart representing the correction procedure. 
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Min = 0.0473 Max = 0,4641 

LAvg= 0,2485 AAvg= 0,2413 

File : , NRay : 100, NView : 100 

Fig. 4.2 First reconstructed image of SI at energy 60 keV. 



Min = 0.0569 Max = 0.4392 

LAvg= 0.2535 AAvg= 0.2415 
File :,NRay : 100, NView : 100 


I'ig. 4.3 Second reconstructed image of SI at energy 60 keV. 


r\ ^ 
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Min = 0.0437 Max = 0.4442 
LAvg = 0.2540 AAvg = 0.2434 
File : . NRav : 100. NView : 100 


!'ig. 4.4 Third reconstructed image of S 1 at energy 60 keV. 



Min = 0.0561 Max = 0.4229 
LAvg = 0,2494 AAvg = 0.2406 
File : , NRay : 100, NView : 100 

Fig.4 . 5 Monchromatic reconstructed image of S 1 at energy 60 keV . 


Correction of beam hardening effect 


Fig. 4.2 is the first reconstructed image of SI obtained from the projection data 
f(p) in which third circle appears light. The area outside the circles, which is green in the 
original image (Fig. 3.1) appears very light. The second reconstructed image (Fig.4.3) is 
obtained from the data m -f(p) + f {p) in which all the three circles appear clearly. The 
third reconstructed image (Fig. 4.4) is obtained from the projection data 
/(7?) + ./'(p)’ matches well with the monochromatic reconstructed image (Fig. 
4.5). The maximum pixel value in the beam hardening reconstmcted image (Fig. 3.2) is 
0.4961, which becomes 0.4442 after the third reconstruction (Fig. 4.4). 
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i Min =-0.0018 Max = 0,9397 
! LAvg = 0.5257 AAvg = 0.4983 
i File : , NRay ; 100, NView : 100 


4.6 F'irst reconstructed image of S2 at energy 0.3 MeV. 



Min = 0.0038 Max= 0.9223 
LAvg = 0.5242 AAvg =0.5016 
File :,NRay : 100, NView ; 100 


}• ig. 4 . 7 Third reconstructed image of S2 at energy 0.3 MeV . 


Correction of beam hardening effect 



: Min = —0.0174 Msx = 0.9120 

I LAvg= 0.5207 AAvg= 0.4974 
! File :,N Ray : 100, NView : 100 


I'ig. 4.8 Moiuicliromalic reconstructed image of S2 at energy 0.3 MeV. 

F'ig. 4.6 i.s the first reconstructed image of S2 obtained from the projection data 
fip). Tliird jvcun.Niructcd image (Fig. 4.7) is obtained from the projection data 
m fi 'p) t ,/'(/>). which i.s matching well with the monochromatic reconstructed image 
(i’ig. 4.8 1. 'file inasinimn pixel value in the beam hardening reconstructed image (Fig. 
3.4) i.s 1 ..3 1 , which becomes 0.9223 after third reconstruction (Fig. 4.7). 
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77 Max = 0.4661 
:310 AAvg= 0.3485 
ly : 100, NView : 100 


rucled image of S3 at energy 60 keV 
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Min = 0.0428 Max = 0.4549 
lAvg = 0,2655 AAvg = 0.3556 
Fite : , NRay : 100, NView : 100 


F'ig. 4. ! I reconstructed image of S3 at energy 60 keV. 

l-'tg. 4.'^> i-< the llrst reconstructed image of S3 obtained from the projection data 
fQh in whiclt t!tc star appears dark blue. The second reconstructed image (Fig.4.10) is 
obtained Inun the data m ■“/(/>) + /(/?) in which star area appears as light blue colour. 
Second rcciniNtiuetcd image (Fig. 4.10) is matching well with the monochromatic 
rectm.stnicicd image (i’ig. 4.11). The maximum pixel value in the beam hardening 
rccon.stntcted image (Fig. .3.6) is 0.6775, which becomes 0.4598 after the second 
reeon.stnietion (Fig, 4.10). 
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Chapter 5 

Results and Discussion 

5J Error iiieasiirenients 

Avcruj'.c pixel cmnx aiul root mean square errors are calculated for projection 
data ax ucll as tta* rcconstructiwi data for all the three specimens. Tables 5.1 to 5.4 give 
the errtn’ estimates for Kl, Tables 5.5 to 5.8 give the error estimates for S2. Tables 5.9 to 
5,12 pjvc the error c},ttmatc,s for STResuUs are quite encouraging. Error values are 
dccrcasinp after each tleratiwi. 
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I'licrgN' 

in kc\' ) 

m-p I 

1 m-f{p) 1 

m-[m-f(p) + f(p)] 1 


0.107.1 

0.0150 

0.0035 

60 

0.0100 

0.0050 

0.0011 

5."' 

0.0.t64 

0.0065 

0.0016 

41 

0.1656 

0.0317 

0.0077 

1 af'k- 

, 1 A\crap,c pixel errors between the projection data for SI . 


h'nciy,> 

m p ] 

j fn~f{p) 

m-[m-f(p) + f{p)] 1 

( in KtA' 1 

i ! 1. I 

u i,'4,' 

omob 

\ 

0.0086 

66 

0.024.1 

! 0.0070 

i ; 

0.0029 


0.0410 

1 0.0086 

0.0038 

41 

; 0.201 1 

! 0,0430 

( 

1 

j 

0.0182 


i al'tk* 5 - Hnat mean squart* errors between the projection data for SI. 
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i 

: j 

i I-'nergy i 

1 - 

Bcjun hardening 

Monochromatic 

After applying 

(inke\’) ; 

i * 



Correction 

101' 

(t.Sdlb 

0.0247 

0.0276 

0(1 

J 

(!.lf4!! 1 

0.0340 

0.0345 

5.1 

nt)5!4 1 

0.0405 

0.0411 

i 41 

f 

1 

t!,115d i 

5 

i 

, w . . 1 

0.0558 

0.0593 

i 


i'ai>U* >. ^ Av'cr.i;!c pi\c! errors between the original and various reconstructed data for 
SI, 


j 


Monocimunatic 

i 

i 

i 

After applying 

Correction 

1 

1 (inkeV'i j 



t 1 1 ’ 

f . * f - # *w 

J 0.0462 

0.0471 

i «’ * 

0,0045 

j 0.0600 

0.0599 

1 

1 5 ^ i 

, f f 1 

i 1 


1 0.0703 

0.0705 

1 41 

1 

1 

‘ 1 

1 

0 

; 0,0977 

1 

\ 

0.1009 


Tiiblc 5 4 Ko»tt ineast square errors between the original and various reconstructed data 

for S! 
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Energy 

( in MeV ) 

— — 

1 fn-p 1 

1 m-f(p) 1 

m-[m-f(p) + f(p)] 1 

0.1 

0.9908 

0.0440 

0.0121 

0.2 

0.1307 

0.0021 

0.0007 

0.3 

0.3383 

0.0057 

0.0017 


'Fable 5.5 Average pixel errors between the projection data for S2. 



Table 5.6 Root mean square errors between the projection data for S2. 
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Energy 

(in MeV) 

Beam hardening 

Monochromatic 

After applying 

Correction 

0.1 

0.6048 

0.2222 

0.2244 

0.2 

0.1581 

0.1008 

0.1009 

0.3 

0.2477 

0.0789 

0.0792 


'Fable 5.7 Average pixel errors between the original and various reconstructed data for 
S2. 


j 

1 

F'ncrgy 

(in MeV) 

Beam hardening 

Monochromatic 

After applying 

Correction 

0.1 

0.6970 

0.3714 

0.3721 

0.2 

0.2050 

0.1705 

0.1706 

0.3 ^ 

0.2658 

0.1340 

0.1342 


Table 5.8 Root mean square errors between the original and various reconstructed data 
for S2. 
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Energy 

( in keV ) 

1 - p j 

i 

1 

1 ^-f{p) 1 

m-[m-f{p) + f{p)] 

100 

0.4945 

0.0277 

0.0175 

84 

0.4235 

0.0236 

0.0150 

60 

0.2343 

0.0138 

0.0091 

52 

0.0157 

0.0027 

1 

0.0017 

41 ! 

0.4787 

0.0391 

0.0238 


I'ablc 5.9 Average pixel errors between the projection data for S3. 


lOicrp) 


m-f{p) 1 

+ /(/?)] 

( in keV ) 




100 

0.5082 

0.0379 

0.0257 

84 

0.4350 

0.0324 

0.0220 

60 

0.2399 

0.0193 

0.0134 

52 

0.0178 

0.0033 

0.0028 

41 

0.4971 

0.0521 

0.0346 


'Fable 5.10 Root mean square errors between the projection data for S3 . 


Results arid Discussions 


Energy 

( in keV) 

Beam hardening 

Monochromatic 

After applying 

Correction 

100 

o.wn 

0.0309 

0.0383 

84 

0.2735 

0.0399 

1 

0.0452 

60 

0.1834 

0.0637 

0.0645 


0.1003 

0.0923 

0.0930 

41 

0.3301 

0.1569 

0.1681 


Table 5.11 Average pixel errors between the original and various reconstructed data for 
S3. 


I- net g} 

Beam hardening 

Monochromatic 

After applying 

(in keV) 



Correction 

1 00 

0.33 14 

0.0553 

0.0581 

84 

0.2936 

0.0702 

0.0706 

60 

0.2063 

0.1101 

0.1089 

52 

0.1638 

0.1577 

0.1585 

41 

0.3650 

0.2656 

0.2756 


Tabic 5.12 Root mean square errors between the original and various reconstructed data 
for H3. 
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5.2 Simulation results 

Due le> tfie beam hardening effect it may be possible that the cross section under 
consideration may he subject to change. But such situation is did not occur in any of the 
three specimens simulated. 

Present simuUition resulLs also show that second and third order polynomials are 
suflicient for correcting the beam hardening effects. After applying the correction the 
rccon.struction results shows that they are not much dependent on the initial guess. 

As for as the choice of the energy of the reconstruction is considered, from the 
error estimates and reconstructed images it is observed, that most of the times the energy 
tor which the root mctin square error between the polychromatic and monochromatic 
projection daut is least, gives the best results. The observed results shows that 2 to 3 
reconstructions are enough for conecting the beam hardening effect, because no 
significant inijir.o ement in the reconstruction was subsequently found. 
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Chapter 6 

Conclusions 


III all the three tested spccimcn.s it is observed that the polychromatic projection 
data is ([iruiTitaiiw!;, different from the monochromatic projection data, but not 
qualitatively. More .'dniulaliun is to be done with different assumed cross sections. 

I'lie reconstructed result.s reported arc quite good in view of the error estimates. 
After applying the concctinn the polychromatic reconstructed images match well with 
the monocliromatic reconstructed images. 

In the pre.sent work the number of rays and the number of views are taken as 100 
in all tfjc three tested .specimen.s. More .simulation has to be done with different number of 
ray.s and views in order to see their effect. 
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